Whole-heart electromechanical simulations using Latent Neural Ordinary Differential Equations

Cardiac digital twins provide a physics and physiology informed framework to deliver personalized medicine. However, high-fidelity multi-scale cardiac models remain a barrier to adoption due to their extensive computational costs. Artificial Intelligence-based methods can make the creation of fast and accurate whole-heart digital twins feasible. We use Latent Neural Ordinary Differential Equations (LNODEs) to learn the pressure-volume dynamics of a heart failure patient. Our surrogate model is trained from 400 simulations while accounting for 43 parameters describing cell-to-organ cardiac electromechanics and cardiovascular hemodynamics. LNODEs provide a compact representation of the 3D-0D model in a latent space by means of an Artificial Neural Network that retains only 3 hidden layers with 13 neurons per layer and allows for numerical simulations of cardiac function on a single processor. We employ LNODEs to perform global sensitivity analysis and parameter estimation with uncertainty quantification in 3 hours of computations, still on a single processor.

Cardiac digital twins integrate physiological and pathological patientspecific data to monitor, analyze and forecast patient disease progression and outcomes.High-fidelity multi-scale and anatomically accurate models are available but require extensive high-performance computing resources to run, which limit their clinical translation 1 .Over the past years, these mathematical models evolved from an electromechanical description of the human ventricular activity in idealized shapes 2,3 and realistic geometries [4][5][6][7] , while also addressing diseased conditions [8][9][10][11][12] , to whole-heart function [13][14][15][16][17][18] .Nevertheless, running many electromechanical simulations still entail high computational costs, hindering the development and application of cardiac digital twins.The use of Machine Learning tools, such as Gaussian Processes Emulators 19 and Artificial Neural Networks (ANNs) 20,21 , allows to create efficient surrogate models that can be employed in many-query applications 22 , such as sensitivity analysis and parameter inference [23][24][25] .In the framework of digital twinning and personalized medicine, bridging the chasm between the need for a supercomputer [26][27][28][29] and performing accurate numerical simulations on a standard computer [30][31][32][33] would have a tremendous impact on the clinical adoption of computational cardiology.
In this work, we develop a Scientific Machine Learning method to build a comprehensive surrogate model involving both cardiac and cardiovascular function.Specifically, we train a system of Latent Neural Ordinary Differential Equations (LNODEs) 20,34,35 that learns the pressurevolume transients of a heart failure patient while varying 43 model parameters that describe cardiac electrophysiology, active and passive mechanics, and cardiovascular fluid dynamics, by employing 400 3D-0D closed-loop electromechanical training simulations.We design a suitable loss function that is minimized during the tuning process of the ANN parameters, which entails small relative errors of LNODEs, i.e., from 2% to 6%, when the number of training samples is small compared to the dimensionality of the parameter space and the explored model variability.These LNODEs enable four-chamber heart numerical simulations on a standard computer by encoding pressure-volume dynamics while spanning electro-mechano-fluid model parameters throughout the cardiovascular system.Furthermore, they can be easily trained on a single central processing unit (CPU).
We use the trained LNODEs to perform global sensitivity analysis (GSA) and robust parameter estimation with uncertainty quantification (UQ) 23,24 .For the former, we observe how model parameters impact the variability of scalar quantities of interest (QoIs) retrieved from the pressurevolume time traces, by considering both first-order and high-order interactions via Sobol indices 36 .For the latter, we combine two Bayesian statistics methods, i.e., Maximum a Posteriori (MAP) estimation and Hamiltonian Monte Carlo (HMC) 34,37,38 , where we exploit efficient matrixfree adjoint-based methods, automatic differentiation and vectorization 34 .In particular, we design several test cases where we calibrate tens of model parameters by matching the pressure and volume time traces, that are timedependent QoIs, coming from 5 unseen 3D-0D numerical simulations for the trained ANN.GSA and parameter estimation with UQ can be carried out in 3 hours of computations by using a single core standard laptop.

Results
We display the whole computational pipeline in Fig. 1.
• Top-left: we use a database of N sims = 405 electromechanical simulations generated by a personalized anatomy four-chamber heart model from a heart failure patient (see Supplementary Material 1), where we vary N P ¼ 43 parameters that describe cell, tissue, whole-heart and cardiovascular system material properties and boundary conditions.For all the numerical simulations, we run 5 heartbeats in sinus rhythm and we perform our analysis on the pressure and volume transients of the last cardiac cycle.We refer to Supplementary Material 2 for all the details about the four-chamber physics-based mathematical model and the numerical settings of these simulations.All the information regarding model parameters can be found in Supplementary Material 3. • Bottom-left: we employ N train, valid = 400 simulations to tune the LNODEs hyperparameters.This surrogate model learns the atrial and ventricular pressure-volume temporal dynamics of the last cardiac cycle only, while receiving time and model parameters as inputs.We perform K-fold cross validation with K = 10 for the training-validation splitting.We detail the whole optimization process to get the final values of the LNODEs hyperparameters in Supplementary Material 4.
We evaluate the accuracy of the trained LNODEs on a testing dataset consisting of the remaining N test = 5 numerical simulations.• Bottom-right: we employ the trained LNODEs to perform GSA.
• Top-right: we estimate model parameters with UQ on N test = 5 numerical simulations by means of the trained LNODEs.
Learning atrial and ventricular pressure-volume loops Automatic hyperparameters tuning with K-fold cross validation leads to an optimal ANN architecture comprising 3 hidden layers and 13 neurons per hidden layer.The optimal number of states is set to N z = 8, i.e., no latent variables are selected.This is motivated by the trade-off between the size of the training set N train, valid with respect to the number of parameters N P , i.e., a thrifty system of LNODEs with no additional hidden variables z latent (t) is selected to avoid overfitting.More details regarding LNODEs training and hyperparameters tuning are given in Supplementary Material 4.
In Table 1, we report the Normalized Root Mean Square Error (NRMSE) and R 2 coefficients associated with the LA, LV, RA, and RV pressure-volume time traces provided by LNODEs.These values are obtained by considering a test set comprised of N test = 5 electromechanical simulations.The accuracy obtained by our surrogate model in reproducing the cardiac outputs is high, manifesting testing errors that approximately range from 2% to 6% for all time-dependent QoIs.The good match between models M 3D-0D and M ANN is also confirmed by Fig. 2, where atrial and ventricular pressurevolume traces present a good overlap on the whole testing set.

Global sensitivity analysis
Figure 3 shows the total-effect Sobol indices.We consider a parameter to be relevant if the associated Sobol indices are greater than Fig. 1 | Sketch of the computational pipeline.We perform several 3D-0D closedloop four-chamber heart electromechanical simulations.We build an accurate and efficient ANN-based surrogate model of the whole cardiovascular function by means of LNODEs.We carry out GSA to understand how each model parameter influences different QoIs extracted from the simulated pressure-volume loops.We robustly estimate many model parameters from time-dependent QoIs.Fully personalized 3D-0D numerical simulations can be performed after parameter calibration with UQ. 10 −1 for at least one QoI.We notice that, as expected from a physiological point of view, some model parameters are compartmentalized, i.e., cell-to-organ level values coming from a certain compartment of the cardiocirculatory system mostly explain the variability of QoIs that are specific to that region.Indeed, some parameters of the CRN-Land model, such as atrial calcium/troponin complex when 50% of crossbridges are blocked perm CRNÀLand , have an important role in determining atrial behavior.Similar considerations occur for the ventricular part of the heart, where the most important parameters are related to the ToRORd-Land model.Nevertheless, it is important to notice the interplay between some ventricular parameters of the ToRORd-Land model at the cellular scale, such as ventricular steadystate duty ratio dr ToRORd-Land , ventricular calcium/troponin complex when 50% of crossbridges are blocked perm ToRORdÀLand 50 and ventricular reference Ca 2+ sensitivity ca ToRORdÀLand 50 and the atrial function.This is a particularly interesting insight into cardiac physiology that can be clearly unraveled using this type of comprehensive sensitivity analysis.We highlight that, as expected, some model parameters, such as atrioventricular delay AV delay , systemic resistance R sys and pulmonary resistance R pulm strongly affect all QoIs, whereas others, such as the pericardial coefficient k peri , as well as aorta parameters (length Aol, stiffness k Art ), have a minor role in determining all QoIs.
Finally, we remark that Sobol indices are affected by the amplitude of the ranges in which the parameters are varied.In particular, the wider the range associated with a parameter, the greater the associated Sobol indices will be, as the parameter in question potentially generates greater variability in the QoI.Therefore, we stress that our results are valid for the specific ranges we used.

Robust parameter estimation
In the context of parameter calibration, a preliminary GSA allows to determine the identifiability of model parameters according to the provided QoIs.Based on the results obtained from global sensitivity analysis, we design 4 in silico test cases to show the robustness and flexibility of our parameter calibration process, which is driven by a combined use of MAP estimation and HMC starting from time-dependent QoIs.In Table 2, we report the observed pressure-volume time traces and estimated model parameters for each test case.In T LV and T ventricles , we estimate model parameters related to the ventricular and cardiovascular function starting from time-dependent QoIs localized in the ventricles.In T atria , we calibrate model parameters over the whole cardiac function and cardiocirculatory network by only considering atrial observations.Finally, we challenge our surrogate model by taking all cardiac pressures and volumes over time and by estimating 11 model parameters.
We perform parameter estimation with UQ on N test = 5 electromechanical simulations that are unseen by the trained LNODEs.Figure 4 shows some two-dimensional views of the posterior distribution for each test case and for all N test numerical simulations.We notice that the true parameter values are contained inside the 95% credibility regions.Moreover, by using Bayesian statistics we are able to capture relationships among model parameters.In particular, in Fig. 4 we consider different pairs of model parameters for each test case and numerical simulation to maximize the number of interactions.For instance, pulmonary resistance scaling factor R pulm and ventricular steady-state duty ratio dr ToRORd-Land are positively correlated with systemic resistance scaling factor R sys and ventricular reference Ca 2+ sensitivity ca ToRORdÀLand 50 , respectively, while ventricular steady-state duty ratio dr ToRORd-Land and ventricular reference Ca 2+ sensitivity ca ToRORdÀLand 50 are negatively correlated with fast endocardial layer scaling factor k FEC and ventricular calcium/troponin complex when 50% of crossbridges are blocked perm ToRORdÀLand 50 , respectively.We notice that, in some cases, cell-based atrial and ventricular parameters may be correlated, as it happens for atrial Ca 2+ -troponin cooperativity TRPN CRNÀLand n and ventricular calcium/troponin complex when 50% of crossbridges are blocked perm ToRORdÀLand 50 , while in most situations, such as with ventricular steady-state duty ratio dr ToRORd-Land and atrial Ca 2+ -troponin cooperativity TRPN CRNÀLand n , there is no interaction.We also remark that this kind of relationships may be unraveled among different physical problems.For instance, this occurs between cardiovascular hemodynamics (systemic resistance scaling factor R sys ) and the ventricular cell tension model (ventricular steady-state duty ratio dr ToRORd-Land ).The aforementioned interactions among model parameters can be quite interesting and surprising from a physiological perspective, especially when they involve very different cardiovascular compartments of this complex multiscale and multiphysics mathematical model.For the sake of completeness, in Table 3, we report the identified parameter values of ventricular steady-state duty ratio dr ToRORd-Land , systemic resistance scaling factor R sys and pulmonary resistance scaling factor R pulm for all test cases, with respect to the first testing simulation.We see that the true values of the parameters are always contained inside the interval defined by median plus/minus interquartile range.We refer to Supplementary Material 6 for the tables containing similar results and comparisons for all test cases (T LV , T ventricles , T atria and T all ) with all the relevant model parameters over the N test electromechanical simulations.

Discussion
In this work, we propose a surrogate model based on LNODEs to learn the pressure-volume temporal dynamics of 3D-0D closed-loop four-chamber heart electromechanical simulations 23 .The geometry is retrieved from a heart failure patient and some functional aspects have been incorporated in the electromechanical simulations.Specifically, we get QRS duration from 12-lead electrocardiograms, ventricle reference tension for active contraction and heart rate, in order to achieve pressure-volume loops that are consistent with measured peak pressure and pressure transient duration 25 .In particular, starting from 400 numerical simulations, we create an anatomy-specific surrogate model by leveraging LNODEs.These are defined by a lightweight feedforward fully-connected ANN containing 3 hidden layers and 13 neurons per layer.LNODEs retain the variability of 43 model parameters that describe electrophysiology, active and passive mechanics, and hemodynamics, both at the cell level and organ scale, and covering a wide range of pressure and volume values (see Figures in Supplementary Material 2).Indeed, LNODEs allows to capture complex dynamics with a small number of tunable parameters.This paradigm, opposed to large machine learning models that work in the overparameterization regime, has proven to be very effective and robust, showing great generalization properties.Some other examples are given by Latent Dynamics Networks 39 or Liquid Neural Networks 40 , which are built on top of LNODEs to account for complex space-time processes exhibiting abrupt changes by using very simple architectures and small latent spaces.
The generation of such a comprehensive training dataset poses an incredible technological challenge itself in the scientific community 25 , and in recent years different surrogate models of cardiac electromechanics based on emulators have been proposed in the literature to provide fast and accurate evaluations based on computationally expensive physics-based mathematical models 19,25,[41][42][43] .These emulators are built on a collection of pre-computed numerical simulations obtained by sampling the parameter space, similarly to what has been done in this work.However, they only fit a static map between model parameters and pointwise QoIs extracted from the numerical simulations.On the other hand, LNODEs present a higher representational power, because they encode time dependent electromechanical simulations instead of pointwise QoIs, while also requiring a smaller amount of data to reach a prescribed accuracy 23 .This is why this paper provides a comprehensive surrogate model embracing cardiac and cardiovascular function.
The choice of 400 numerical simulations for the training of LNODEs allows to get consistently low validation and testing errors, in the order of 2% to 6%, even in areas of the parameter space that are sparsely covered by the training samples (see Figure in Supplementary Material 3).The error remains within these bounds even if we increase the dimension of the testing : VE calcium/troponin complex when 50% of crossbridges are blocked, dr ToRORd-Land : VE steady-state duty ratio, ca ToRORdÀLand 50 : VE reference Ca 2+ sensitivity, CV ventricles : VE conduction velocity in the fiber direction, k FEC : fast endocardial layer scaling factor, R sys : systemic resistance scaling factor, R pulm : pulmonary resistance scaling factor, AT: atrial, VE: ventricular.
set while decreasing the one of the training set (see Supplementary Material 9).In addition, LNODEs provide better generalization properties compared to Gaussian processes emulators, especially for ventricular function (see Supplementary Material 8).
LNODEs require a small amount of computational resources and enable several applications of interest in a very fast and accurate manner.Indeed, as reported in Table 4, running the training phase of the ANN along with GSA and robust parameter estimation on a single core standard laptop just requires 13 h of computations.We remark that this time can be reduced with a multi-core implementation.On the other hand, employing the 3D-0D model M 3D-0D for the same computational pipeline would entail very significant costs.The overall speed-up with the surrogate model M ANN is equal to 1718x.Furthermore, after the parameter calibration process, a whole-heart electromechanical simulation can run via high-performance computing with the estimated model parameters, showing all relevant space-time fields, such as transmembrane potential, active contraction, displacement and stresses.This would provide relevant insights across the whole high-fidelity cardiac model.
It is interesting to note that several model parameters related to electrophysiology, mechanics, and hemodynamics at the cell-to-organ scale have a significant impact on the pressure-volume loops.These model parameters can be inferred from the pressure-volume relationships using Bayesian parameter estimation with UQ, i.e., their true values are contained within the 95% credibility regions of the posterior distribution.In addition, bayesian statistics provides important insights by capturing crosscorrelations between model parameters.This occurs even between different cardiovascular compartments, such as the systemic circulation and ventricular electromechanics.Our approach could be applied in a clinical setting by using clinically measured rather than computational pressurevolume loops to infer protein and cell-to-organ function directly from clinical data.
Even though this paper focuses on a single anatomy, it represents an important milestone towards the construction of emulators incorporating geometric variability.Indeed, the presented approach can be extended to cover patient variability, by incorporating statistical shape modeling 44,45 or other ANN-based methods, such as Universal Solution Manifold Networks 46 or generative deep learning techniques based on Signed  Distance Fields 47,48 , within LNODEs, to encode different geometrical parameterizations.Furthermore, multiple pathological conditions and diagnoses can be taken into account by providing specific one-hot vectors as additional inputs to the ANN 47 .In this manner, we would run the numerical simulations with the biophysically detailed and anatomically accurate mathematical model just once and we would train an ANN that generalizes on multiple patients, while effectively capturing multi-physics and multiscale knowledge.Moreover, even though patient-specific pressure-volume loops have not been considered in this work, we aim at adding this information as part of our computational pipeline for parameter calibration with UQ.The proposed method paves the way to extensions incorporating different anatomies and pathological conditions, which would potentially allow for a universal whole-heart simulator that might be readily deployed in clinical practice for fast and reliable personalized parameter calibration based on patient-specific data.

Ethics statement
The clinical data used in this study were collected as part of a clinical trial (REC reference: The ANN receives N z state variables, N P scalar parameters, and two periodic inputs.Indeed, even though LNODEs are just trained on the last cardiac cycle, the cosine and sine terms account for the heartbeat period T HB and the atrioventricular delay AV delay of whole-heart electromechanical simulations (see Supplementary Material 2 for further details).On the other hand, the vector of physics-based model parameters θ, which involves cardiac electromechanics and cardiovascular hemodynamics, is not related to the time variable, as is the case for T HB and AV delay , and is given directly as input neurons to the ANN.We stress that, differently from ref. 23, the initial reduced state vector z 0 contains different sets of initial conditions for pressures, volumes and latent variables.Pressure and volume initial values are determined by model M 3D-0D .Following 49 , latent variables are initialized to zero and these initial conditions act as additional tunable parameters along with the weights and biases of the ANN.
The loss function that we minimize during the ANN optimization process reads: w for the ANN.It comprises the normalized mean square error between ANN pressure-volume predictions z physical (t) and observations zphysical ðtÞ, as well as a weak penalization of the physical state vector time derivatives, maximum and minimum values for t ∈ [T − T HB , T].Indeed, given the small ratio between the dimensionality of the training dataset and the number of parameters θ of model M 3D-0D , we notice that these three additional terms reduce the generalization errors of the ANN.The penultimate weakly enforced condition on z latent (t) favors a periodic solution for all the hidden latent variables.The last term of the loss function prescribes the L 2 regularization of the ANN weights and ι is one of the automatically tuned LNODEs hyperparameters (see Supplementary Material 4).

Global sensitivity analysis
We employ the Saltelli's method to perform a variance-based sensitivity analysis 50 .We compute both first-order Sobol indices and totaleffect Sobol indices for each combination of quantity of interest and model parameter 36 .These two indices define how much varying a single parameter affects a specific QoI and how higher-order interactions among model parameters influences the model outputs, respectively.All mathematical details regarding the computation of

50 ,
atrial Ca 2+ -troponin cooperativity TRPN CRNÀLand n and atrial reference Ca 2+ sensitivity ca CRNÀLand 50 , or of the Guccione model, such as atrial stiffness in the transverse plane b atria t

Fig. 2 |
Fig. 2 | Pressure and volume time transients obtained with model M 3D-0D (dashed lines), compared to those obtained with model M ANN (solid lines), on the testing samples (N test = 5).Light blue: LA, orange: LV, blue: RA, green: RV.

Fig. 4 |
Fig. 4 | Two-dimensional views of the posterior distribution estimated by means of HMC for each test case (rows) and N test = 5 electromechanical simulations (columns).Different colors are associated to T LV , T ventricles , T atria , T all .

Table 1 |
Testing errors and R 2 coefficients on the timedependent outputs of the trained LNODEs system

Table 2 |
Summary of the 4 in silico test cases for parameter calibration LA , p RA , p LV , p RV , V LA , V RA , V LV , V RV dr ToRORd-Land , perm ToRORdÀLand

Table 3 |
True value and median with interquartile range (between brackets) associated to the estimated values of ventricular steady-state duty ratio dr ToRORd-Land , systemic resistance scaling factor R sys and pulmonary resistance scaling factor R pulm during HMC for the first numerical simulation of the testing set 14/WM/1069) approved by the West Midlands-Coventry & Warwickshire Research Ethics Committee.and ventricular active tension or passive stiffness, and resistances of the systemic and pulmonary circulation.The reduced state vector zðtÞ 2 R N z contains the time-dependent pressure and volume variables of the left atrium (LA), right atrium (RA), left ventricle (LV) and right ventricle (RV), as well as additional latent variables without a direct physical interpretation, that is zðtÞ ¼ ½z physical ðtÞ; z latent ðtÞ ¼ ½p LA ðtÞ; p LV ðtÞ; p RA ðtÞ; p RV ðtÞ; V LA ðtÞ; V LV ðtÞ; V RA ðtÞ; V RV ðtÞ; z latent ðtÞ T .
:ð1Þwhere z 0 is the vector of initial conditions.The ANN, with weights and biases encoded in w 2 R N w , is defined by ANN : R N z þ2þN P !R N z .Vector θ 2 Θ & R N P defines the model M 3D-0D parameters.Some examples of θ could be conductances of different ionic channels, myocardial conductivity, atrial

Table 4 |
Summary of the approximated computational times to perform GSA and parameter estimation with UQ. 3D-0D closed-loop model M 3D-0D (top) and LNODEs M ANN (bottom)